*Main data pathway globals, set directory to clean data folder
global results "S:\Project\DemoSos2\common\felles\JR_RG\DrVA\ResultsRev\"
global cleandata "S:\Project\DemoSos2\common\felles\JR_RG\DrVA\CleanData\"
cd "$cleandata"

*Make new folder and subfolders 
capture mkdir "$cleandata"
capture mkdir "$results"

cap log close
clear all
set more off
set matsize 11000
set maxvar 120000


use  clean_patientlevel_file_rev3.dta, clear
keep if yr_str_exog_swap >= 2005 & yr_str_exog_swap<=2014
keep if str_exogGPIDnew!=. 


* largest mobility group 
a2group, individual(str_exogGPIDnew) unit(str_exogGPIDprev) groupvar(pair)
bys pair: ge size_pair = _N
tab pair
drop if size_pair < 10 ///largest mobility group has 99%

keep if str_exog_age >=55


ge death5 = mortality_5year
ge death2 = mortality_2year
for num 25 35 45 55 65: ge death2_X = mortality_2year if str_exog_age>=X
for num 1 3 4 5: ge deathX_55 = mortality_Xyear if str_exog_age>=55
for var heart cancer ext resp mental accident suicide homicide: ge X2y_55 = X_2y if str_exog_age>=55

ge str_exogGPIDprev_yearswp = yr_str_exog_swap*str_exogGPIDprev

label var death2_55 	"2-years mortality for 55+"
label var bweight_2y 	"Birthweight"
label var heart2y_55 	"Mortality: heart conditions"
label var cancer2y_55 	"Mortality: cancer"
label var ext2y_55 		"Mortality: external conditions"
label var resp2y_55 	"Mortality: respiratory conditions"
label var mental2y_55 	"Mortality: mental health conditions"
label var accident2y_55	"Mortality: accident conditions"

replace total_base = 25 if total_base >= 25 & total_base!=.

*** get a sample of doc that dont retire (stop before retirement or stop working as GP)
ge noretirement = str_exog_drage_prev < 60 | min_age_prev== max_age_prev

* patient counts
ge d = death2_55!=.
bys str_exogGPIDnew: egen count_pat = sum(d)
drop d


/* estimate the FE */
local var "death2_55"
* base
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var'_base= mean(fe2_`var') 
ge n`var'_base = e(N) if e(sample)	
bys str_exogGPIDnew: egen s = sum(fe2_`var')
ge d = 1 if e(sample)
bys str_exogGPIDnew: egen n = sum(d) 
ge lofe_`var'_base = (s-fe2_`var')/(n-1)
drop _fe_`var' fe2_`var' res_`var'  n d s

* mun FE
reghdfe `var' i.yr_str_exog_swap i.str_exog_age male i.str_exog_muni, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var'_mun = mean(fe2_`var') 
ge n`var'_mun = e(N) if e(sample)	
bys str_exogGPIDnew: egen s = sum(fe2_`var')
ge d = 1 if e(sample)
bys str_exogGPIDnew: egen n = sum(d) 
ge lofe_`var'_mun = (s-fe2_`var')/(n-1)
drop _fe_`var' fe2_`var' res_`var'  n d s


* old GPxyear of SWAP FE
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev_yearswp) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var'k_t = mean(fe2_`var') 
ge n`var'k_t = e(N) if e(sample)	
bys str_exogGPIDnew: egen s = sum(fe2_`var')
ge d = 1 if e(sample)
bys str_exogGPIDnew: egen n = sum(d) 
ge lofe_`var'k_t = (s-fe2_`var')/(n-1)
drop _fe_`var' fe2_`var' res_`var'  n d s



* include baseline characteristics
local baseline "hs0 NORborn antbarn tot_income dwoverfor bigcity birthorder i.total_base"
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age `baseline', a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var'_pre = mean(fe2_`var') 
ge n`var'_pre = e(N) if e(sample)	
bys str_exogGPIDnew: egen s = sum(fe2_`var')
ge d = 1 if e(sample)
bys str_exogGPIDnew: egen n = sum(d) 
ge lofe_`var'_pre = (s-fe2_`var')/(n-1)
drop _fe_`var' fe2_`var' res_`var'  n d s

*Panel B
* drop from estimation if less than X count_pat
foreach num in 25 50 100 { 
	reghdfe `var' i.yr_str_exog_swap male i.str_exog_age if count_pat >= `num', a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
	ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
	bys str_exogGPIDnew: egen fe_`var'_c`num' = mean(fe2_`var') 
	ge n`var'_c`num' = e(N) if e(sample)	
	drop _fe_`var' fe2_`var' res_`var'
}


* Panel C
local var "death1_55"
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var' = mean(fe2_`var') 
ge n`var' = e(N) if e(sample)	
drop _fe_`var' fe2_`var' res_`var'

local var "death3_55"
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var' = mean(fe2_`var') 
ge n`var' = e(N) if e(sample)	
drop _fe_`var' fe2_`var' res_`var'

local var "death4_55"
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var' = mean(fe2_`var') 
ge n`var' = e(N) if e(sample)	
drop _fe_`var' fe2_`var' res_`var'

local var "death5_55"
reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var' = mean(fe2_`var') 
ge n`var' = e(N) if e(sample)	
drop _fe_`var' fe2_`var' res_`var'

* Panel D 
local var "death2_65"
reghdfe `var'  i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var') 
ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
bys str_exogGPIDnew: egen fe_`var'_65 = mean(fe2_`var') 
ge n`var'_65 = e(N) if e(sample)	
drop _fe_`var' fe2_`var' res_`var'


	
bys str_exogGPIDnew: gen count_doc = _n
local var "death2_55"
twoway 	(kdensity fe_`var'_base if count_doc==1, lcolor(red) lpattern(dash) lwidth(thick) graphregion(color(white)) bgcolor(white) xtitle(Value Added)) ///
		(kdensity fe_`var'_docy if count_doc==1, lcolor(blue) lpattern(solid) lwidth(thick)) ///
		(kdensity fe_`var'_pre if count_doc==1, lcolor(black) lpattern(solid) lwidth(thick)), legend(order(1 "Baseline VA" 2 "No gender/age FE" 3 "Extra Controls VA")) 
		graph export  "$results\figure_a6.pdf",  replace	
	
*---------------------------------------------------------*
* OUTPUT INFO
*---------------------------------------------------------*

/* TABLES */
global fe "fe_death2_55_base fe_death2_55_mun fe_death2_55k_t fe_death2_55_pre fe_death1_55 fe_death3_55 fe_death4_55 fe_death5_55 "///
	"fe_death2_55_c25 fe_death2_55_c50 fe_death2_55_c100 fe_death2_65"

label var fe_death2_55_base 	"baseline"
label var fe_death2_55_mun 		"municipality FE"	
label var fe_death2_55k_t 		"pre-swap GP-year FE"
label var fe_death2_55_pre 		"pre-swap chronic conditions"
label var fe_death2_55_c25 		"25+ patients"
label var fe_death2_55_c50 		"50+ patients"
label var fe_death2_55_c100		"100+ patients"
label var fe_death1_55 			"1-years mortality for 55+"
label var fe_death3_55 			"3-years mortality for 55+"
label var fe_death4_55 			"4-years mortality for 55+"
label var fe_death5_55 			"5-years mortality for 55+"
label var fe_death2_65			"2-years mortality for 65+"

eststo s1: estpost tabstat $fe, s(N p10 p25 p50 p75 p90 mean sd) columns(statistics)  
esttab s1 using "$results/fe_unadjusted_rob.tex", replace fragment ///
	cells("p10(fmt(3)) p25(fmt(3)) p50(fmt(3)) p75(fmt(3)) p90(fmt(3)) mean(fmt(3)) sd(fmt(3))") noobs nomtitle ///
	varlabels (`e(labels)')  varwidth(20) booktabs label ///
	collabels("N" "P10" "P25" "P50" "P75" "P90" "Mean" "SD")

keep  str_exogGPIDnew fe_* count_pat lopenr lo*
save "$results/fe_addon_rob_lo2.dta", replace 
keep str_exogGPIDnew fe_* count_pat
duplicates drop
save "$results/fe_addon_rob2.dta", replace

eststo s1: estpost tabstat $fe, s(N) columns(statistics)  
esttab s1 using "$results/fe_unadjusted_rob_N.tex", replace fragment ///
	cells("count(fmt(0))") noobs nomtitle ///
	varlabels (`e(labels)')  varwidth(20) booktabs label 

	
	
use "$results/fe_addon_rob2.dta", clear
for var fe_death1_55 fe_death2_55_base fe_death3_55 fe_death4_55 fe_death5_55: egen st_X= std(X)

foreach var in st_fe_death1_55 st_fe_death3_55 st_fe_death4_55 st_fe_death5_55 {
	eststo `var': reg `var'  st_fe_death2_55_base, robust
}
esttab st_fe_death1_55 st_fe_death3_55 st_fe_death4_55 st_fe_death5_55 using "$results/table_a13.tex",  replace keep(st_fe_death2_55_base) ///
	b(4) se(4) star(* 0.1 ** 0.05 *** 0.01) stats(emp N, fmt(%9.3g)) 


	
*-------------------------------------------------------------------------
* new mobility group restriction - based on patients age (Panel E, Table A12)
*-------------------------------------------------------------------------
use  clean_patientlevel_file_rev3.dta, clear
keep if yr_str_exog_swap >= 2005 & yr_str_exog_swap<=2014
keep if str_exogGPIDnew!=. 

keep if str_exog_age >=55
* largest mobility group 
a2group, individual(str_exogGPIDnew) unit(str_exogGPIDprev) groupvar(pair)
bys pair: ge size_pair = _N
tab pair
drop if size_pair < 10 ///largest mobility group has 99%


sort pair
egen gr = group(pair)
tab gr
drop gr

ge death2 = mortality_2year
for num 55: ge death2_X = mortality_2year if str_exog_age>=X

label var death2_55 	"2-years mortality for 55+"

/* estimate the FE */
global index "death2_55"
foreach var in   $index {
	reghdfe `var' i.yr_str_exog_swap male i.str_exog_age, a(_fe_`var'=str_exogGPIDnew str_exogGPIDprev) res(res_`var')  

	ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 	
	bys str_exogGPIDnew: egen fe_`var' = mean(fe2_`var') 
	ge n`var'_base = e(N) if e(sample)	
	drop _fe_`var' fe2_`var' res_`var'
	
}
keep str_exogGPIDnew fe_death2_55* 
duplicates drop

*---------------------------------------------------------*
* OUTPUT INFO
*---------------------------------------------------------*
label var fe_death2_55 	"2-years mortality for 55+: strict mobility group"

eststo s1: estpost tabstat fe_death2_55, s(count p10 p25 p50 p75 p90 mean sd) columns(statistics)  
esttab s1 using "$results/table_a12.tex", append fragment ///
	cells("N(fmt(3)) p10(fmt(3)) p25(fmt(3)) p50(fmt(3)) p75(fmt(3)) p90(fmt(3)) mean(fmt(3)) sd(fmt(3))") noobs nomtitle ///
	varlabels (`e(labels)')  varwidth(20) booktabs label ///
	collabels("N" "P10" "P25" "P50" "P75" "P90" "Mean" "SD")

	


	
